Beyond the locality approximation in the standard diffusion Monte Carlo method 
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We present a way to include non local potentials in the standard Diffusion Monte Carlo method 
without using the locality approximation. We define a stochastic projection based on a fixed node 
effective Hamiltonian, whose lowest energy is an upper bound of the true ground state energy, even 
in the presence of non local operators in the Hamiltonian. The variational property of the resulting 
algorithm provides a stable diffusion process, even in the case of divergent non local potentials, like 
the hard-core pseudopotentials. It turns out that the modification required to improve the standard 
Diffusion Monte Carlo algorithm is simple. 



PACS numbers: 02.70.Ss, 31.10.+Z, 31.25.-v 



Diffusion Monte Carlo (DMC) is one of the most suc- 
cessful methods to compute the ground state properties 
of quantum systems. Although the fixed node (FN) ap- 
proximation is needed to cure the infamous sign problem 
for fermions, the accuracy of the DMC framework has 
yielded many benchmark results Q. However, when the 
DMC method is applied to "ab initio" realistic Hamilto- 
nians, its computational cost scales cx Z , where Z is 
the atomic number |2j. Therefore, the use of pseudopo- 
tentials is necessary to make those calculations feasible. 

Since the pseudopotentials are usually non local, the 
"locality approximation" is made besides the FN, by re- 
placing the true Hamiltonian H with an effective one 
H eS , which reads 0: 



K + Vi 



J dx'{x'\V non loc \x)^ T (x') 



(1) 



where K is the kinetic operator, V\ oc is the local poten- 
tial, and the last term in Eq. ^is the non local potential 
localized by means of the trial wave function vPy. The 
projection is then realized by iteratively applying the op- 
erator G = exp(— t(H cS — E c g)) to 9t in order to filter 
out its high energy components. The localized potential 
enters in the branching part (birth and death process) of 
the algorithm, while the usual FN constraint is employed 
to limit the diffusion process within the nodal pockets of 
^r, and avoid the fcrmionic sign problem. Thus E e f{ is 
the FN ground state energy of H cS , computed during the 
sampling of the mixed distribution \]/ e g^/y: 



(y oS \H cii \y T ) (* eff |If|* T ) 



(*cff|*T) 



(*cff|*T) 



= E M a- (2) 



Ema is the mixed average of H, and the above identity 
holds because H eS ^ T /^ T = H^ T /^ T - Since * e ff is the 
FN ground state of H cS , which differs from H, Ema is no 
longer equal to the variational FN energy of H , defined 
as: 



E F N = <*cff|#|*cff)/(*cff|*cff>- 



(3) 



Therefore, in contrast with the case of local Hamilto- 
nians, Ema calculated with the locality approximation 



does not in general give an upper bound to the ground 
state energy of H (variational principle) . 

In a previous work0, we introduced the Lattice Reg- 
ularized Diffusion Monte Carlo algorithm (LRDMC), 
which provides an upper bound for the true ground state 
energy and allows estimate Efn, even in the case of non 
local potentials. In this paper we propose an extension 
of the standard DMC framework that gives the same re- 
sults as the LRDMC method, after a proper modification 
of the DMC propagator. 

We start by considering the importance sampling 
Green function 



G(x <— x, t) 



x '\ e -r(H-E T , 



\*), (4) 



where Et is an energy offset, r the time step, and x 
a vector of particle coordinates. In the diffusion Monte 
Carlo method, G(x' <— x, r) is iteratively applied to ^1., 
in order to sample stochastically the mixed distribution 
$(x,t) = ^t(x)^(x ) t), *if(x,t) converging to the lowest 
possible state in energy. To rewrite G(x' x, r) (Eq.0J 
in a practical way, it is necessary to resort to the Trotter 
break up, which is exact in the limit of r — » 0. Here we 
split the Hamiltonian into local and non local operators, 
and we end up with the following expression for the Green 
function: 



G(x' ^x,t)~ J dx" T x , !X „(t) G DMC (x" «- x,t), (5) 
where Gdmc{ x ' x i T ) is the usual DMC propagator^, 



(2 7 rr) i 



exp 



(x' — X — Tv{x)Y 

2t 



(6) 

and T x i jX (t) is the matrix containing the non local po- 
tential, 

£zV) tf\ e -Tv m l0C | x) ^ §x _ tVx , (7) 

In the above Eqs. N is the total number of parti- 
cles, v{x) = Vln|* T (a;)| the drift velocity, E l £ c (x) = 



2 



(K + V\ oc )^>t{x) I^t{x) the contribution to the local 
energy coming from the local operators, and V x ' lX = 
% T T %) ( x '\ V non ioc\x). The final form of Gdmc has been 
obtained by further splitting the Hamiltonian into the 
kinetic and potential part, while the exponential of the 
non local potential in T has been linearized up to order 

T. 

If the case of pseudopotentials, the number of non- 
zero matrix elements V x ', x will be finite, once a quadra- 
ture rule with a discrete mesh of points is applied to 
evaluate the projection over the angular components 
of the pseudopotential^, Therefore, the process in 
G(x' <— x, t) driven by T x / )X (r) can be calculated using a 
heat bath algorithm, since T x > tX {r)/ ^2 X » T x h >x (t) can be 
seen as a transition probability, and it can be computed 
a priori for all possible new coordinates x' . We notice 
that the matrix elements T x ' tX (r) are easily evaluated in 
a standard DMC algorithm, since 14',x are already com- 
puted to calculate the localized pseudopotential in Eq. ^ 



/ dx'{x'\V non i oc \x)V T (x') 



'Vx> 



(8) 



At variance with the locality approximation, V x ', x con- 
tribute now to move the particles, according to the tran- 
sition matrix T (T-moves). 

An important limitation of this idea is given by the 
sign problem. Indeed both ^fey and (a;'|V^ lon i oc |^) can 
change sign, which should be included in the weights, but 
this yields averages with exponentially increasing noise. 
A solution is to apply the FN approximation not only to 
Gdmc but also to T, which becomes: 



T™(t) = 6 x , x ,-tV x I x) 



(9) 



where we defined x = 1/2 (14', a; ± |Vx',x|)- In prac- 
tice, we keep only those matrix elements which give a 
positive T x t >x [r). Moreover, we add to the diagonal po- 
tential the so called "sign flip term" , i.e. the sum over 
the discarded matrix elements V x , . Therefore, the local 
potential becomes 



V° s (x) = V loc (x) + J2 V ^- 



(10) 



This is equivalent to work with a new effective FN Hamil- 
tonian 



H°% = K + V cS (x) 

rreff 

x' ,x 



(11) 



(x'\Vnon loc\x) if V x >,x < 0. 



In contrast to the effective Hamiltonian of the locality 
approximation written in Eq. ^ the ground state energy 
E eS (= E M a) of the above H eS is an upper bound for the 
ground state energy of the true H. As shown in Ref. 6 



for the Lattice Green function Monte Carlo, this varia- 
tional property is due to the sign flip term (positive con- 
tribution) added to the local potential, and the T-moves 
driven by the off diagonal matrix elements V x , x . Instead, 
in the locality approximation also V~, is summed in the 
diagonal part (Eq. |5J) , and this leads to an attractive po- 
tential, which cannot provide a variational property for 
Em a- Moreover, we found that the negative divergences 
of the fully localized potential on the nodes of ^>t are re- 
sponsible in some case (e.g. see Fig. ^) for numerical in- 
stabilities in the locality approximation, which disappear 
once H in Eq.^Jis used together with the T^^-moves. 
Indeed, whenever V~, x is large, it pushes the walker away 
from the attractive regions of the localized potential, and 
protects the sampling from divergences in the weights. 
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FIG. 1: Energies for the Carbon pseudoatom with r = 
0.08-ff -1 at the given DMC generation, is an antisym- 
metrized geminal pow er (AGP) wave function with a 3-body 
Jastrow factor |l3 . Ha|. We report results for the locality ap- 
proximation (H an with a — 1 and 7 = 0) and the algorithm 
with T-moves (a — 0,7 = 0). 

Once a T FN -move is generated according to the 
transition probability T^ n x (t)/ £ x „ T£, n x (t), the walker 
should acquire the weight wt{x,t) — Yl x " ^x'^x( T ) ^ ue 
to the normalization of the T FN matrix. This weight can 
be recast as an exponential form valid up to order r, 



Wt 



exp 



/ j x,x 



(12) 



Thus the overall weight w(x, r) of G(x' <— x, t) will be 

w(x, t) = w D mc w T = exp [-t(E l (x) ~ E T )} , (13) 

where wdmc is the weight of Gdmc for the effec- 
tive Hamiltonian (V\ oc replaced by V cfi ), and El(x) — 
H eS ^ T /^ T = H^t/^t is the local energy. Notice that 
a non-symmetric branching factor has been included in 
Gdmc (Eq. EJ). When we use the exponential form in 
Eq. El and consequently the weight in Eq. ^| the time 
step error is usually smaller than that obtained with the 
linear form. This can be understood in the limit of per- 
fect importance sampling. Indeed, if is close to the 



3 



ground state of H, the weight in Eq. ^] is almost con- 
stant, since the variance of El(x) is small, and the time 
step bias is reduced. 

The proposed DMC scheme for fixed node Hamiltoni- 
ans with non local potentials is the following: (i) perform 
a diffusion-drift move according to G<us(x' x, t) = 
exp [—(a;' — x — tv(x)) 2 /2t] /(2ttt)~ as is done in the 
standard DMC algorithm, and accept or reject this move 
according to the probability 



,{x') 



G diS (x' <- x,t)^(x) 



(14) 



(ii) weight the walker with the factor 
exp [— t(El{x') — Et)}; (hi) displace the walker a second 
time, with a T-move selected according to the transition 
probability p(x" - x>, r) = T™,(r)/ J2 y T™(r), com- 
puted a priori for all possible new x" . The branching 
process will be the same as in the usual DMC algorithm. 
In practice, only the T-move is the new step, which 
is performed after weighting the walker 0- Although 
we perform an acceptance/rejection step (Eq. | 14H . 
which has been shown to reduce the time step error |8| 
in GdmCi the algorithm does not satisfy exactly the 
detailed balance except in the limit of r — > 0, due to the 
break up of G into Gdmc and T FN (Eq.EJl, and the use 
of a non symmetric branching factor in Eq. El 

In order to estimate the variational FN energy Efn 
(Eq. |2J), and study the quality of the locality ap- 
proximation, we introduce a more general effective 
Hamiltonian[3| H a ^, 



h«;2 = K + v loc (x) + (i + y)J2 v ^ 

x' 

+a(l + 7)ZX« 



(15) 



H%1 = -7 (x'\V non loc \x) if V x >, x > 

Hpl = (1 - a(l + 7 )) {x'\V non loc \x) if V x ,, x < 0, 

where < a < 1 and < 7 < 1/a — 1 are two exter- 
nal parameters. In order to sample the Green function 
G(x' <— x,t) for H ar< , it is sufficient to modify the ma- 
trix T x i -x (t), which becomes 



T Q ' 7 = 

x' ,x 



T 7 Ft 

/ X ,x 



r (l-a(l + 7 )) V- 



if x 

if K 

if Vx 



x' 

X 1 .X 



< 0. 



(16) 



The ground state E(a, 7) of H aa is equal to Ema{&-, 7) 
(Eq.0), since H a ^^ T /^T = H^t/^t by construction. 
The Hamiltonian in Eq. 1111 is recovered with a = and 
7 = 0, while the Hamiltonian of the locality approxima- 
tion (Eq.^l is obtained with a = 1 and 7 = 0. Therefore, 
H an can interpolate between these two extremes, but the 
variational principle for Ema(<1)7) is not guaranteed as 
soon as a ^ 0, since the attractive term a(l+7) J2x< x 



is added to the diagonal potential. However by means of 
H an one can estimate the value of Epiy(a,j) (Eq. 01, 
which is variational for every a and 7, since it is the ex- 
pectation value of the true H on the ground state of H a, t. 
Indeed H = H a ^ - (1 + i)d 1 H a -~< , and the Hellmann- 
Feynman theorem leads to the relation 



E FN (a,j) = E{a,i) - (1 +7) d 1 E(a,i). 



(17) 



One can show 9] that, for a given value of a, the lowest 
Efn(o>,j) is obtained for 7 = 0. Therefore, in order to 
find the best variational estimate of the ground state of 
H , it is enough to calculate the expression in Eq. 1171 with 
7 = 0. In this way one can check which a provides the 
best variational state for H . The derivative djE(a, 0) 
can be computed with either finite differences or corre- 
lated sampling. In both cases, one should keep in mind 
that 7< 1/a — 1, to guarantee the positivity of the T aa 
matrix fEa. I16fl . and so calculating Ep^{a,Q) becomes 
harder as a gets closer to 1. 

Here we present the application of the method to 
the Si and C pseudoatoms. We computed Ema{cx,0) 
and E FN (a,0) for a = 0,0.5,0.9, and the DMC en- 
ergy with the locality approximation, which corresponds 
to -Ema(1,0). With an aim to quantify the locality er- 
ror, and the correction provided by the effective Hamil- 
tonian _ff Q ' 7 , we used three different trial wave func- 
tions (with no Jastrow, a two-body, and a three-body 
(electron-electron- ion) Jastrow factor respectively), shar- 
ing the same dctcrminantal part, and hence the same 
nodes. In this way, the FN error can be separated from 
the effect of the locality approximation, which causes a 
dependence of the DMC energy on the shape of ^>t- 

For the Si atom we used an s — p norm-conserving 
Hartree-Fock (HF) pseudopotential, which is soft and 
has been generated using the Vanderbilt construction [Io|. 
The determinantal part of is a HF wave function with 
a 6s6p/l s|p_ Gaussian basis set. The 2-body Jastrow is 
from Ref. |llj, while the 3-body Jastrow factor is from 
Ref. ^3 Both of the Jastrow factors have been opti- 
mized using an energy minimization procedure ^3 ■ The 
results are reported in Fig. [5] The variational Ep]y(a, 0) 
improves going from a = to a = 0.9, i.e. approach- 
ing the locality approximation. It means that at least 
for this soft pseudopotential the locality approximation 
(a = 1,7 = 0) gives a ground state which is a good vari- 
ational wave function for H. Notice however that the 
standard DMC energies Ema{1,0) have a sizable local- 
ity error, while E F n with a = 0.9 depends only slightly 
on the shape of the trial wave function. A similar re- 
sult was obtained with the LRDMC method for the same 
pseudoatomQ. 

For the C atom we chose to work with an SBK 
pseudopotential 14], which is extremely hard, since it 
diverges like 1/r 2 in the s channel, and 1/r in its lo- 
cal component. The Slater part of is an antisym- 
metrized geminal power (AGP) wave function[ll 13 
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FIG. 2: (Color online) Ema(q,0) (green, dashed line) and 
Efn{ol,0) (blue, dotted line) energies for the Silicon pseu- 
doatom with different values of the effective Hamiltonian pa- 
rameter q. The DMC energies with locality approximation 
(red, solid line), corresponding to Ema(1,0), are reported in 
all panels for reference. We used three different ^t's, which 
have the same determinantal part. A more accurate cor- 
responds to a smaller difference between its variational energy 
(Evmc) and the Bma(0, 0) energy, reported in abscissa. 




FIG. 3: (Color online) The same as Fig. [5] but for the Carbon 
pseudoatom. 

with a 2s2p Gaussian basis set, optimized in the presence 
of the 3-body Jastrow factor by minimizing its variational 
energy 0| . The determinantal part has been kept fixed 
in the other two ^t's, which differ only by their Jastrow 
factors. The results are plotted in Fig. [3] Here the local- 
ity approximation is very poor, as it leads to non varia- 
tional Em a- The spikes in Fig. ^ coming from regions 
of the configuration space where the effective potential is 
attractive, are surely responsible of the non variational 
results. Surprisingly, without Jastrow, which has a 
higher energy, leads to much more stable DMC simula- 
tions. The locality approximation, which relies on the 
quality of the shape of in the core, performs poorly 
with this hard-core pseudopotential, since it is difficult 



to find the optimal shape of in the core region, due 
to the divergence of the non local pseudopotential. In- 
deed EpN(a, 0) is higher for a = 0.9, being the worst for 
the 3-body Jastrow factor. On the other hand, the best 
variational Efn(&,0) is obtained for a = 0, irrespective 
of the form of the Jastrow factor. 

To summarize, we have described a scheme to treat non 
local potentials within the standard DMC method. We 
have extended the DMC formalism to handle a generic 
Hamiltonian with discrete off-diagonal matrix elements 
and the fixed node approximation. Only a simple modi- 
fication of the standard algorithm is required to include 
the T-moves generated according to the non local po- 
tentials. By using an effective Hamiltonian approach, we 
showed that it is possible to have stable simulations, even 
in the case of divergent hard-core pseudopotentials, and 
obtain variational results. A similar effective Hamilto- 
nian has been successfully used in the LRDMC method. 
The difference is in the kinetic part, which is discretized 
in the lattice regularized approach. The LRDMC and 
the DMC methods have the same efficiency for small Z, 
although it is possible to have a gain in the LRDMC 
efficiency by an ad hoc choice of the kinetic parameters, 
particularly for heavier elements. We conclude, by noting 
that the same Green function presented here can be used 
in the Reptation Quantum Monte Carlo method 

This work was supported by the NSF grant DMR- 
0404853. We thank S. Sorella, C. Umrigar, D. M. Ceper- 
ley, S. Chiesa, and J. Kim for useful discussions. 



[1] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Ra- 

jagopal, Rev. Mod. Phys. 73, 33 (2001). 
[2] D. M. Ceperley, J. Stat. Phys. 43, 815 (1986); A. Ma, 

N. D. Drummond, M. D. Towler, and R. J. Needs, Phys. 

Rev. E 71, 066704 (2005). 
[3] L. Mitas, E. L. Shirley and D. M. Ceperley, J. Chem. 

Phys. 95, 3467 (1991). 
[4] M. Casula, C. Filippi, and S. Sorella, Phys. Rev. Lett. 

95, 100201 (2005). 
[5] S. Fahy, X. W. Wang and Steven G. Louie, Phys. Rev. B 

42, 3503 (1990). 
[6] D. F. B. ten Haaf, H. J. M. van Bemmel, J. M. J. van 

Leeuwen, W. van Saarloos, and D. M. Ceperley, Phys. 

Rev. B 51, 13039 (1995). 
[7] Since the normalization of T™ x f Eg. 1121 depends on x, 

the best choice is to weight the walker before the T-move. 
[8] P. J. Reynolds, et al, J. Chem. Phys. 77, 5593 (1982); 

C. J. Umrigar, et al, J. Chem. Phys. 99, 2865 (1993). 
[9] S. Sorella, cond-mat/0201388 
[10] D. Vanderbilt, Phys. Rev. B 32, 8412 (1985). 
[11] C. Filippi and C. J. Umrigar, J. Chem. Phys. 105, 213 

(1996). 

[12] M. Casula, C. Attaccalite, and S. Sorella, J. Chem. Phys. 

121, 7710 (2004). 
[13] S. Sorella, Phys. Rev. B 71, 241103 (2005). 
[14] W. J. Stevens, H. Basch, and M. Krauss, J. Chem. Phys. 



81, 6026 (1984). 

M. Casula, and S. Sorella, J. Chem. Phys. 119, 6500 
(2003). 



[16] S. Baroni, and S. Moroni, Phys. Rev. Lett. 82, 4745 
(1999). 



